#include <stdio.h>
#include <stdlib.h>
#include <complex.h>
#include <math.h>
#include <fftw3.h>

#define N 44100

#include "gnuplot_i.h"


typedef double trame[N];

/******************gnuplot**********************/
gnuplot_ctrl * plot_init(char * style){
  gnuplot_ctrl * h = gnuplot_init();
  gnuplot_setstyle(h, style);
  return h;
}

void plot_temporal(gnuplot_ctrl *h, trame x_axis,trame y_axis, int n, char * nom){
  gnuplot_resetplot(h);
  gnuplot_plot_xy(h, x_axis, y_axis, n, nom);
  sleep(80);
}

/******************sinusoide**********************/
double * generer_sinusoide(int n){
  double * tab = malloc(n * sizeof(double));
  int i;
  for(i=0; i < n; i++)
    tab[i] = 0.5 * sin(2.0 * 3.141592654 * 440.0 * (double)i/44100.0);

   for(i=0; i < n; i++)
    tab[i] += 0.5 * sin(2.0 * 3.141592654 * 1200.0 * (double)i/44100.0);

  return tab;
}

/******************dft**********************/
void dft(double s[N], complex S[N]){
  const double pi = 4 * atan(1.0);
  complex c2 = 0;
  int m, n;
  for(m=0; m < N; m++){
    c2 = 0;
    for(n=0; n < N; n++)
      c2 += s[n] * cexp(-I * (2*pi/N) * n*m);
    S[m] = c2;
  }
      
}

/**************fft************************/



void fft(trame s, fftw_complex data_out[N])
{
  static fftw_plan plan;
 fftw_complex data_in[N];

 plan = fftw_plan_dft_1d(N,data_in,data_out, 
			 FFTW_FORWARD, 
			 FFTW_ESTIMATE);

  int i;
  for(i = 0; i < N; i++)
    data_in[i] = s[i];
  fftw_execute(plan);
}

int main(){
  double * sinusoide = generer_sinusoide(N);
  complex cplx[N];
  //dft(sinusoide, cplx);
  fft(sinusoide, cplx);
   int i;
  double ordonee[N];
  for(i=0; i < N; i++)
  ordonee[i] = cabs(cplx[i]);//module du complexe

 
    

  double abscisse[N];
  for(i=0; i < N; i++)
    abscisse[i] = i ;//*44100.0/N;

  gnuplot_ctrl * h = plot_init("lines");
  plot_temporal(h, abscisse, ordonee, N/4, "dft");
  
  free(sinusoide);
  return 0;
}
